Meta-analysis of the Interoceptive Accuracy Scale (IAS)

Study 1

library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Data Preparation

Sample 1

• Murphy (2020)
  • https://osf.io/3m5nh

df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

Sample 2

• Gaggero (2021)
  • https://osf.io/5x9sg

load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  mutate(across(starts_with("IAS "), as.numeric)) |>
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)

Sample 3

• Campos (2022) - Study 1
  • https://osf.io/j6ef3

df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Sample 4

• Todd (2022)
  • https://osf.io/ms354/

df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Sample 5

• Arslanova (2022)
  • https://osf.io/mp3cy/
  • note: Final.xlsx includes information about which participants passed the attention checks

df5_attention <- readxl::read_excel("C:/Users/asf25/Documents/studies/InteroceptionIAS/data/Arslanova2022/Participants2.xlsx", 
    sheet = "FINAL") |>
  filter(ActorHR == 1) |>
  select(Subj.ID, Gender, Age) |>
  mutate(Gender = as.numeric(Gender), Age = as.numeric(Age)) |>
  filter(!is.na(Age)) |> #error on the documentation of this participant - attention check failed but not noted
  mutate(Gender = case_when(
    Gender== 0 ~ "Male",
    Gender== 1 ~ "Female" # based on paper reporting 65 women
  ))

New names:
• `Why No?` -> `Why No?...4`
• `Why No?` -> `Why No?...6`
• `` -> `...12`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`

df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |>
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  rename(Subj.ID = "Participant.Public.ID") |>
  select(-starts_with("Participant")) 

df5 <- df5 |>
  filter(Subj.ID %in% df5_attention$Subj.ID) |>
  select(-Subj.ID) |>
  mutate(across(everything(), as.numeric))

Sample 6

• Brand (2022)
  • https://osf.io/xwz6g/

df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Sample 7

• Brand (2023)
  • https://osf.io/3f2h6/

df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER, -Sample) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
  )

Sample 8

• Lin (2023)
  • https://osf.io/3eztd
  • Note: No tickle because same Chinese character

df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH)

Sample 9

• VonMohr (2023) - Study 3
  • https://osf.io/7p9u5/

df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Sample 10

• Makowski
  • Note: First sample has some missing items
  • Note: Analog scales
  • note df10c = primals data

df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    typeSample = Sample
  )

df10c <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionPrimals/refs/heads/main/data/data_participants.csv") |>
  select("Participant" = "participant_id", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |>
  filter(!is.na(Heart)) # participant is outlier and had total IAS scores below 0.4 

pi_vars <- c("PI_Enticing", "PI_Alive", "PI_Safe", "PI_Good", "PI_Changing", 
             "PI_Hierarchical", "PI_Understandable", "PI_AboutMe", "PI_Abundant", 
             "PI_Acceptable", "PI_Beautiful", "PI_Cooperative", "PI_Funny", 
             "PI_Harmless", "PI_Improvable", "PI_Intentional", "PI_Interconnected", 
             "PI_Interesting", "PI_Just", "PI_Meaningful", "PI_NeedsMe", 
             "PI_Pleasurable", "PI_Progressing", "PI_Regenerative", "PI_Stable", 
             "PI_WorthExploring")


correlation::correlation(
  df10c[, "IAS_Total", drop = FALSE],  
  data2 = select(df10c, all_of(pi_vars)), 
) |> 
  # Formatting correlation results
  mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
  
  # Releveling factors for visualization
  mutate(Parameter2 = fct_relevel(Parameter2, rev(pi_vars)),
         Parameter1 = fct_relevel(Parameter1, "IAS_Total")) |>
  
  # Plotting the correlation matrix
  ggplot(aes(x = Parameter1, y = Parameter2)) +
  geom_tile(aes(fill = r), color = "white") +  
  geom_text(aes(label = val), size = 3) +  
  labs(title = "Correlation Matrix") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", limit = c(-1, 1), midpoint = 0, guide = guide_colourbar(ticks = FALSE)) +  
  theme_minimal() +
  theme(
    legend.title = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(hjust = 1)
  )

Sample 11

• Poreiro et al.,
  • Study 1: https://osf.io/49wbv
  • Study 2: unpublished

df11a <- haven::read_sav("../data/Giulia/Interoceptive Attention ESM.sav") |> 
  select(AGE, GENDER, contains("IAS_ACC"), -IAS_ACC_20, starts_with(c("MAIA", "TAS20", "DEP", "ANX")),  -ANX_21, -TAS20_21) |> #anx_21 , tas20_21 are attention checks
  mutate_all(as.numeric) |>
  rename(
    Age = AGE, Gender = GENDER,
    Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4,
    Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7, Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9,
    Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
    Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18,
    Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_21, Itch = IAS_ACC_22
  ) |>
   mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender %in% c(3, 5) ~ "Other",
    Gender == 4 ~ "NA"
  )) |> ## compute dimensions for the TAS
  mutate(TAS_DIF = TAS20_1 + TAS20_3 + TAS20_6 + TAS20_7 + TAS20_9 + TAS20_13 + TAS20_14,
         TAS_DDF = TAS20_2 + (1-TAS20_4) + TAS20_11 + TAS20_12 + TAS20_17,
         TAS_EOT = (1-TAS20_5) + TAS20_8 + (1-TAS20_10) + TAS20_15 + TAS20_16 + (1-TAS20_18) + (1-TAS20_19) + TAS20_20)


df11b <- haven::read_sav("../data/Giulia/Sleep and Exteroceptive Interoceptive Sensitivity.sav") |> # 165 participants
  filter(Total_Attention_Fail == 0) |> # 32 participants removed based on failing one or more checks
  select(AGE, GENDER, contains("IAS_ACC"), ) |> 
  mutate_all(as.numeric) |> 
  rename( Age = AGE, Gender = GENDER, 
          Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4, Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7,
          Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9, Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
          Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18, Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_20, 
          Itch = IAS_ACC_21 )  |> 
  mutate(Gender = case_when(
    Gender == 1 ~ "Male", 
    Gender == 2 ~"Female", 
    Gender %in% c(3, 5) ~ "Other", 
    Gender == 4 ~ "NA", )) |>
  na.omit()


# Save all
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b, df10c=df10c, df11a=df11a, df11b=df11b)
save(data, file = "../data/data.RData")

Participants

• Sample 1a: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
• Sample 1b: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
• Sample 2: Data from Gaggero’s(2020) study, downloaded from OSF, included 814 participants (Mean age = 24.9, SD = 5.3, range: [18, 58], 0.2% missing; Gender: 60.3% women, 39.4% men, 0.25% non-binary).
• Sample 3: Data from Campos’s(2022) study, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
• Sample 4: Data from Todd’s(2022) study, downloaded from OSF, included 802 participants ().
• Sample 5: Data from Arslanova (2022) study, downloaded from OSF, included 143 participants (Mean age = 28.5, SD = 7.6, range: [18, 73]; Gender: 45.5% women, 54.5% men, 0.00% non-binary).
• Sample 6: Data from Brand’s(2022) study, downloaded from OSF, included 619 participants (Mean age = 43.9, SD = 14.5, range: [18, 78]; Gender: 78.7% women, 20.2% men, 1.13% non-binary).
• Sample 7a: Data from Brand’s(2023) study, downloaded from OSF, included 522 participants (Mean age = 23.4, SD = 6.7, range: [18, 79]; Gender: 79.5% women, 19.7% men, 0.77% non-binary).
• Sample 7b: Data from Brand’s(2023) study, downloaded from OSF, included 1993 participants (Mean age = 32.0, SD = 12.6, range: [16, 81]; Gender: 77.7% women, 21.7% men, 0.60% non-binary).
• Sample 7c: Data from Brand’s(2023) study, downloaded from OSF, included 808 participants (Mean age = 27.3, SD = 9.3, range: [18, 72], 0.5% missing; Gender: 68.9% women, 30.2% men, 0.87% non-binary).
• Sample 8a: Data from Lin’s(2023) study, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
• Sample 8b: Data from Lin’s(2023) study, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
• Sample 9: Data from VonMohr’s(2023) study 3, downloaded from OSF, included 21843 participants (Mean age = 56.5, SD = 14.4, range: [18, 93], 0.2% missing; Gender: 73.2% women, 25.1% men, 1.55% non-binary, 0.15% missing).
• Sample 10a: Data from [Makowski’s(2023)] study , downloaded from OSF, included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).
• Sample 10b: Data from [Makowski’s(2023)] study , downloaded from OSF, included 836 participants (Mean age = 25.1, SD = 11.3, range: [18, 76], 0.1% missing; Gender: 73.8% women, 22.6% men, 2.87% non-binary, 0.72% missing; Education: Bachelor, 22.85%; Doctorate, 2.15%; High School, 63.52%; Master, 6.22%; missing, 0.24%; Other, 5.02%).
• Sample 10c: Data from [Makowski’s(2023)] study , downloaded from OSF, included 104 participants (Mean age = 21.6, SD = 5.0, range: [18, 50], 1.0% missing; Gender: 76.0% women, 19.2% men, 2.88% non-binary, 1.92% missing).
• Sample 11a: Data from Giulia’s (2024) study , included 107 participants (Mean age = 26.8, SD = 9.2, range: [18, 57]; Gender: 74.8% women, 23.4% men, 1.87% non-binary)
• Sample 11b: Data from [Giulia] study , included 131 participants (Mean age = 30.9, SD = 12.0, range: [18, 60]; Gender: 76.3% women, 22.9% men, 0.76% non-binary)

Total N = 32214.

desc <- function(df){
  summary <- df |> summarise(
      n = n(),
      mean = mean(Age, na.rm = TRUE),
      sd = sd(Age, na.rm = TRUE),
      min = min(Age, na.rm = TRUE),
      max = max(Age, na.rm = TRUE),
      female = sum(Gender == "Female", na.rm = TRUE)
    )
    return(summary)
}

desc_df1a <- desc(df1a)
desc_df1b <- desc(df1b)
desc_df2 <- desc(df2)
desc_df3 <- desc(df3)
#desc_df4 <- desc(df4)
desc_df5 <- desc(df5_attention)
desc_df6 <- desc(df6)
desc_df7a <- desc(df7a)
desc_df7b <- desc(df7b)
desc_df7c <- desc(df7c)
desc_df8a <- desc(df8a)
desc_df8b <- desc(df8b)
desc_df9 <- desc(df9)
desc_df10a <- desc(df10a)
desc_df10b <- desc(df10b)
desc_df10c <- desc(df10c)
desc_df11a <- desc(df11a)
desc_df11b <- desc(df11b)

desc_total <- rbind(desc_df1a, desc_df1b, desc_df2, desc_df3, desc_df5, desc_df6, desc_df7a, desc_df7b, desc_df7c, desc_df8a, desc_df8b, desc_df9, desc_df10a, desc_df10b, desc_df10c, desc_df11a, desc_df11b)

## calculate weighted mean

total_mean = sum(desc_total$mean * desc_total$n) / sum(desc_total$n)
total_sd = sum(desc_total$sd * desc_total$n) / sum(desc_total$n)
perce_female = sum(desc_total$female)/sum(desc_total$n) *100

library(gt)

# APA style ####
gt_apastyle <- function(gt_table, font.size=12) {
  gt_table  |> 
    gt::opt_table_lines(extent = "none") %>%
    gt::tab_options(
      heading.border.bottom.width = 2,
      heading.border.bottom.color = "black",
      heading.border.bottom.style = "solid",
      table.border.top.color = "black",
      table.border.top.style = "solid",
      table.border.top.width = 2,  
      table_body.hlines.color = "white",
      table_body.border.top.color = "black",
      table_body.border.top.style = "solid",
      table_body.border.top.width = 2,
      heading.title.font.size = font.size,
      table.font.size = font.size,
      heading.subtitle.font.size = font.size,
      table_body.border.bottom.color = "black",
      table_body.border.bottom.width = 2,
      table_body.border.bottom.style = "solid",
      column_labels.border.bottom.color = "black",
      column_labels.border.bottom.style = "solid",
      column_labels.border.bottom.width = 1
    ) |> 
      gt::opt_table_font(font = "times")
}

make_age <- function(age) {
  age <- as.numeric(age) 
  
  mean(age, na.rm = TRUE) |> 
    insight::format_value(digits=1) |> 
    paste0(" ± ", 
           insight::format_value(sd(age, na.rm = TRUE), digits=1)) 
} 

# Create the dataframe
table <- data.frame(
   Sample = c('Murphy et al., (2020)', '', '', 'Gaggero et al., (2021)', 'Campos et al., (2022)', 'Todd et al., (2022)', 'Arslanova et al., (2022)', 'Brand et al., (2022)', 'Brand et al., (2023)', '', '', '', 'Lin et al., (2023)', '', '', 'VonMohr et al., (2023)', 'Makowski et al., (2023a)', 'Makowski et al., (2023b)', 'Makowski et al., (2023c)', 'Poreiro et al., (2024)', 'Poreiro et al., unpublished', 'Total'),
  Subsample = c('', 'Sample 1', 'Sample 2', '', '', '', '', '', '', 'Sample 1', 'Sample 2', 'Sample 3', '', 'Sample 1', 'Sample 2', '', '', '', '', '','', ''),   
  Language = c('', 'English', 'English' , 'English and Italian', 'Portuguese', 'English', 'English', 'German', '', 'German', 'German', 'German', '', 'Chinese', 'Chinese', 'English', 'English', 'English', 'English', 'English', 'English',''),
  N = c('', nrow(df1a), nrow(df1b), nrow(df2), nrow(df3), nrow(df4), nrow(df5_attention), nrow(df6),'', nrow(df7a), nrow(df7b), 802,'', nrow(df8a), nrow(df8b), nrow(df9), nrow(df10a), nrow(df10b), nrow(df10c), nrow(df11a), nrow(df11b), 32216),
  Age = c('', make_age(df1a$Age), make_age(df1b$Age), make_age(df2$Age), make_age(df3$Age), "48.6.6 ± 14.1*", make_age(df5_attention$Age), make_age(df6$Age),'', make_age(df7a$Age), make_age(df7b$Age), make_age(df7c$Age),'', make_age(df8a$Age), make_age(df8b$Age), make_age(df9$Age), make_age(df10a$Age), make_age(df10b$Age), make_age(df10c$Age), make_age(df11a$Age), make_age(df11b$Age), '48.6 ± 13.1'),
  Range = c('', '18-69', '18-91', '18-58', '18-72', '18-92*', '18-73', '18-78', '', '18-79', '16-81', '18-72','', '16-60', '20-60', '18-93', '18-73', '17-76','18-50', '18-57', '18-60', '17-93'),
  Female_Percentage = c('', '69.4%', '70.1%', '60.3%', '59.6%', '50%*', '46.8%', '78.7%', '', '79.5%', '77.7%', '68.9%', '', '57.0%', '56.2%', '73.2%', '50.3%', '53.0%', '76%', '74.8%', '75.9%', '71.6'),
  Availability= c('', '', 'https://osf.io/3m5nh', 'https://osf.io/5x9sg', 'https://osf.io/j6ef3', 'https://osf.io/ms354', 'https://osf.io/mp3cy', 'https://osf.io/xwz6g', '', '', 'https://osf.io/3f2h6', '','','', 'https://osf.io/3eztd', 'https://osf.io/7p9u5', 'https://github.com/RealityBending/IllusionGameReliability', 'https://github.com/DominiqueMakowski/PHQ4R', 'https://github.com/RealityBending/InteroceptionPrimals', 'https://osf.io/49wbv', '','')
) 

table_apa <- table |> 
  gt() |>
  cols_align(align = c("right"), columns = "Age") |> 
  cols_label(Age = "Age (Mean  ± SD)", Female_Percentage = "Female %") |> 
  tab_footnote("* Information taken from the sample description of relevant paper rather than recomputed.") |> 
  gt_apastyle() 
  
table_apa
gtsave(table_apa, "figures/table1.tex")

Descriptive

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10a",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10b",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10c <- estimate_density(select(df10c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10c",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens11a <- estimate_density(select(df11a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11a")
dens11b <- estimate_density(select(df11b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11b")


p1 <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b, dens10c, dens11a, dens11b) |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b", "Sample 10c", "Sample 11a", "Sample 11b")) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  scale_color_material_d() +
  scale_linetype_manual(values = c("solid", "dashed"), guide="none") +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        plot.title = element_text(face="bold")) +
  labs(title = "Distribution of responses for all items across various datasets", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free_y", nrow = 5)

ggsave("figures/Figure1.png", p1, width=10, height=10, dpi=200, bg="white")
p1

data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))
data10c <- select(df10c, all_of(dens10c$Parameter))
data11a <- normalize(select(df11a, all_of(dens11a$Parameter)), range = c(1, 5))
data11b <- na.omit(normalize(select(df11b, all_of(dens11b$Parameter)), range = c(1, 5))) # remove NAS
row.names(data11b) <- NULL



data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), mutate(data8b, Tickle = NA), data9,
  mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA),  data10b,  data10c, data11a, data11b
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

All data

make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Parameter1, y = Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}

make_correlation(data_all)

Sample 1a

make_correlation(data1a)

Sample 1b

make_correlation(data1b)

Sample 2

make_correlation(data2)

Sample 3

make_correlation(data3)

Sample 4

make_correlation(data4)

Sample 5

make_correlation(data5)

Sample 6

make_correlation(data6)

Sample 7a

make_correlation(data7a)

Sample 7b

make_correlation(data7b)

Sample 7c

make_correlation(data7c)

Sample 8a

make_correlation(data8a)

Sample 8b

make_correlation(data8b)

Sample 9

make_correlation(data9)

Sample 10a

make_correlation(data10a)

Sample 10b

make_correlation(data10b)

Sample 10c

make_correlation(data10c)

Sample 11a

make_correlation(data11a)

Sample 11b

make_correlation(data11b)

Graph Analysis

Exploratory Graph Analysis (EGA) is a recently developed framework for psychometric assessment, that can be used to estimate the number of dimensions in questionnaire data using network estimation methods and community detection algorithms, and assess the stability of dimensions and items using bootstrapping. See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

All data

uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.364

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.290
 Urinate Defecate 0.265

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.230
  Heart Breathing 0.225
 Hungry   Thirsty 0.217

uva$keep_remove

$keep
[1] "Itch"

$remove
[1] "Tickle"

Sample 1a

uva <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241

uva$keep_remove

NULL

Sample 1b

uva <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"

Sample 3

uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208

uva$keep_remove

$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"

Sample 4

uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.415
 Tickle      Itch 0.351

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.293
 Urinate Defecate 0.282
    Wind     Burp 0.275
  Sneeze    Cough 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i node_j   wto
 Muscles   Pain 0.205

uva$keep_remove

$keep
[1] "Heart" "Itch" 

$remove
[1] "Breathing" "Tickle"   

Sample 5

uva <- EGAnet::UVA(data = data5, cut.off = 0.3)

Warning: Some variables did not belong to a dimension: Blood_Sugar, Affective_touch 

Use caution: These variables have been removed from the TEFI calculation

uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

    node_i node_j   wto
    Tickle   Itch 0.247
    Bruise   Itch 0.246
 Breathing  Vomit 0.207

uva$keep_remove

NULL

Sample 6

uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237

uva$keep_remove

$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  

Sample 7a

uva <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 

Sample 7b

uva <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237

uva$keep_remove

$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 

Sample 7c

uva <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218

uva$keep_remove

$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"

Sample 8a

uva <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213

uva$keep_remove

NULL

Sample 8b

uva <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206

uva$keep_remove

NULL

Sample 9

uva <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24

uva$keep_remove

$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 

Sample 10a

uva <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209

uva$keep_remove

NULL

Sample 10b

uva <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

NULL

Sample 10c

uva <- EGAnet::UVA(data = data10c, cut.off = 0.3)

Warning: Some variables did not belong to a dimension: Heart, Vomit, Blood_Sugar 

Use caution: These variables have been removed from the TEFI calculation

uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

NULL

Sample 11a

uva <- EGAnet::UVA(data = data11a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.330
   Heart Breathing 0.314
    Wind      Burp 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i  node_j   wto
 Muscles    Pain 0.279
  Tickle    Itch 0.277
  Hungry Thirsty 0.255

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
  Vomit Sneeze 0.201

uva$keep_remove

$keep
[1] "Heart"    "Defecate" "Burp"    

$remove
[1] "Breathing" "Urinate"   "Wind"     

Sample 11b

uva <- EGAnet::UVA(data = data11b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Hungry Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j   wto
   Wind    Burp 0.214
   Temp Muscles 0.202

uva$keep_remove

NULL

Networks

All data

make_egas <- function(data) {
  egas <- list()
  for (model in c("glasso")) {  # , "TMFG"
    for (algo in c("walktrap", "louvain")) {
      for (type in c("ega", "ega.fit")) { # "" # "hierega"
        if (type == "ega.fit" & algo == "louvain") next # Too slow
        egas[[paste0(model, "_", algo, "_", type)]] <- EGAnet::bootEGA(
          data = data,
          seed = 123,
          model = model,
          algorithm = algo,
          EGA.type = type,
          type = "resampling",
          plot.itemStability = FALSE,
          verbose = FALSE
        )
      }
    }
  }

  p <- EGAnet::compare.EGA.plots(
    egas$glasso_walktrap_ega, 
    # egas$glasso_walktrap_ega.fit,
    egas$glasso_louvain_ega, 
    # egas$TMFG_louvain_ega,
    # egas$TMFG_walktrap_ega, 
    # egas$TMFG_walktrap_ega.fit,
    labels = c(
      "glasso_walktrap_ega", 
      # "glasso_walktrap_ega.fit",
      "glasso_louvain_ega" 
      # "TMFG_louvain_ega",
      # "TMFG_walktrap_ega", 
      # "TMFG_walktrap_ega.fit"
    ),
    rows = 3,
    plot.all = FALSE
  )$all
  list(egas = egas, p = p)
}


egas_all <- make_egas(data_all)
egas_all$p

Sample 1a

egas1a <- make_egas(data1a)
egas1a$p

Sample 1b

egas1b <- make_egas(data1b)
egas1b$p

Sample 2

egas2 <- make_egas(data2)
egas2$p

Sample 3

egas3 <- make_egas(data3)
egas3$p

Sample 4

egas4 <- make_egas(data4)
egas4$p

Sample 5

egas5 <- make_egas(data5)
egas5$p

Sample 6

egas6 <- make_egas(data6)
egas6$p

Sample 7a

egas7a <- make_egas(data7a)
egas7a$p

Sample 7b

egas7b <- make_egas(data7b)
egas7b$p

Sample 7c

egas7c <- make_egas(data7c)
egas7c$p

Sample 8a

egas8a <- make_egas(data8a)
egas8a$p

Sample 8b

egas8b <- make_egas(data8b)
egas8b$p

Sample 9

egas9 <- make_egas(data9)
egas9$p

Sample 10a

egas10a <- make_egas(data10a)
egas10a$p

Sample 10b

egas10b <- make_egas(data10b)
egas10b$p

Sample 10c

egas10c <- make_egas(data10c)

Warning: Some variables did not belong to a dimension: Heart, Vomit, Blood_Sugar 

Use caution: These variables have been removed from the TEFI calculation

Warning in order(as.numeric(names(dimension_stability))): NAs introduced by
coercion
Warning in order(as.numeric(names(dimension_stability))): NAs introduced by
coercion

egas10c$p

Sample 11a

egas11a <- make_egas(data11a)
egas11a$p

Sample 11b

egas11b <- make_egas(data11b)
egas11b$p

Structure Stability

All data

patchwork::wrap_plots(lapply(egas_all$egas, plot), nrow = 3)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

Sample 1a

patchwork::wrap_plots(lapply(egas1a$egas, plot), nrow = 3)

Sample 1b

patchwork::wrap_plots(lapply(egas1b$egas, plot), nrow = 3)

Sample 2

patchwork::wrap_plots(lapply(egas2$egas, plot), nrow = 3)

Sample 3

patchwork::wrap_plots(lapply(egas3$egas, plot), nrow = 3)

Sample 4

patchwork::wrap_plots(lapply(egas4$egas, plot), nrow = 3)

Sample 5

patchwork::wrap_plots(lapply(egas5$egas, plot), nrow = 3)

Sample 6

patchwork::wrap_plots(lapply(egas6$egas, plot), nrow = 3)

Sample 7a

patchwork::wrap_plots(lapply(egas7a$egas, plot), nrow = 3)

Sample 7b

patchwork::wrap_plots(lapply(egas7b$egas, plot), nrow = 3)

Sample 7c

patchwork::wrap_plots(lapply(egas7c$egas, plot), nrow = 3)

Sample 8a

patchwork::wrap_plots(lapply(egas8a$egas, plot), nrow = 3)

Sample 8b

patchwork::wrap_plots(lapply(egas8b$egas, plot), nrow = 3)

Sample 9

patchwork::wrap_plots(lapply(egas9$egas, plot), nrow = 3)

Sample 10a

patchwork::wrap_plots(lapply(egas10a$egas, plot), nrow = 3)

Sample 10b

patchwork::wrap_plots(lapply(egas10b$egas, plot), nrow = 3)

Sample 10c

patchwork::wrap_plots(lapply(egas10c$egas, plot), nrow = 3)

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_text()`).

Sample 11a

patchwork::wrap_plots(lapply(egas11a$egas, plot), nrow = 3)

Sample 11b

patchwork::wrap_plots(lapply(egas11b$egas, plot), nrow = 3)

Scores

ega <- egas_all$egas$glasso_louvain_ega$EGA
plot(ega)

get_scores <- function(data, ega) {
  EGAnet::net.scores(data, ega)$scores$std.scores |>
    as.data.frame() |>
    setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))
}

ega_scores_1a <- get_scores(data1a, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_1b <- get_scores(data1b, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_2 <- get_scores(data2, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_3 <- get_scores(data3, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_4 <- get_scores(data4, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_5 <- get_scores(data5, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_6 <- get_scores(data6, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_7a <- get_scores(data7a, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_7b <- get_scores(data7b, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_7c <- get_scores(data7c, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_8a <- get_scores(data8a, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_8b <- get_scores(data8b, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_9 <- get_scores(data9, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_10a <- get_scores(data10a, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_10b <- get_scores(data10b, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_10c <- get_scores(data10c, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_11a <- get_scores(data11a, ega)

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

ega_scores_11b <- na.omit(get_scores(data11b, ega)) # remove two rows with NAS

The default 'loading.method' has changed to "revised" in {EGAnet} version >= 2.0.7.

 For the previous default (version <= 2.0.6), use `loading.method = "original"`

row.names(ega_scores_11b) <- NULL

Factor Analysis

How Many Factors

All data

n <- parameters::n_factors(data_all, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 1a

n <- parameters::n_factors(data1a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 1b

n <- parameters::n_factors(data1b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
5	BIC (adjusted)	BIC

Sample 2

n <- parameters::n_factors(data2, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
4	BIC	BIC
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree

Sample 3

n <- parameters::n_factors(data3, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
3	CNG	CNG
3	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE
4	BIC	BIC

Sample 4

n <- parameters::n_factors(data4, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE
5	BIC	BIC

Sample 5

n <- parameters::n_factors(data5, n_max = 12)

Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
1	BIC	BIC
3	CNG	CNG
3	BIC (adjusted)	BIC
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE

Sample 6

n <- parameters::n_factors(data6, n_max = 12)

Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 7a

n <- parameters::n_factors(data7a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 7b

n <- parameters::n_factors(data7b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
5	Bentler	Bentler
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree

Sample 7c

n <- parameters::n_factors(data7c, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Scree (SE)	Scree_SE

Sample 8a

n <- parameters::n_factors(data8a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree
5	BIC	BIC

Sample 8b

n <- parameters::n_factors(data8b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	BIC	BIC
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree
5	Scree (SE)	Scree_SE
5	BIC (adjusted)	BIC

Sample 9

n <- parameters::n_factors(data9, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 10a

n <- parameters::n_factors(data10a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	BIC	BIC
4	BIC (adjusted)	BIC

Sample 10b

n <- parameters::n_factors(data10b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 10c

n <- parameters::n_factors(data10c, n_max = 12)

Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
1	BIC	BIC
3	Bartlett	Barlett
3	CNG	CNG
4	beta	Multiple_regression
5	Optimal coordinates	Scree
5	Scree (SE)	Scree_SE

Sample 11a

n <- parameters::n_factors(data11a, n_max = 12)

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
The estimated weights for the factor scores are probably incorrect.  Try a
different factor score estimation method.

Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully

Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
The estimated weights for the factor scores are probably incorrect.  Try a
different factor score estimation method.

Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
1	BIC	BIC
3	Bentler	Bentler
3	CNG	CNG
3	VSS complexity 2	VSS
4	beta	Multiple_regression
5	Optimal coordinates	Scree

Sample 11b

n <- parameters::n_factors(data11b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

EFA

All data

efa5 <- parameters::factor_analysis(data_all, n = 4, rotation = "oblimin", sort = TRUE)

Loading required namespace: GPArotation

plot(efa5)

display(efa5)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Burp	0.74	0.03	-0.02	-0.07	1.02	0.48
Cough	0.69	-0.02	0.03	0.06	1.02	0.48
Wind	0.67	0.02	-0.02	1.92e-03	1.00	0.54
Sneeze	0.63	-0.05	0.05	0.14	1.13	0.52
Vomit	0.43	0.06	0.04	0.19	1.43	0.63
Temp	0.26	0.25	0.05	0.21	3.01	0.61
Sex_arousal	0.25	0.23	0.02	0.20	2.94	0.67
Taste	0.24	0.20	0.10	0.14	2.99	0.71
Breathing	0.05	0.61	-0.04	0.03	1.03	0.59
Hungry	-0.11	0.54	-5.65e-03	0.20	1.35	0.63
Heart	0.08	0.49	0.03	-0.07	1.11	0.73
Thirsty	-0.07	0.47	-5.43e-03	0.25	1.59	0.64
Pain	0.15	0.40	0.10	0.10	1.55	0.62
Muscles	0.30	0.37	0.07	0.02	2.01	0.59
Blood_Sugar	0.11	0.35	0.21	-0.18	2.42	0.76
Bruise	0.24	0.34	0.22	-0.19	3.26	0.65
Tickle	-0.04	-0.04	0.82	0.04	1.02	0.38
Itch	0.02	5.03e-03	0.75	-0.04	1.01	0.43
Affective_touch	0.07	0.20	0.33	0.11	2.04	0.70
Urinate	0.04	0.10	0.04	0.65	1.06	0.45
Defecate	0.16	0.04	0.05	0.60	1.15	0.47

The 4 latent factors (oblimin rotation) accounted for 41.58% of the total variance of the original data (MR1 = 14.45%, MR3 = 11.76%, MR2 = 8.09%, MR4 = 7.28%).

Sample 1a

efa4 <- parameters::factor_analysis(data1a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR2	MR3	MR1	MR4	Complexity	Uniqueness
Tickle	0.72	0.03	-0.10	0.11	1.09	0.47
Itch	0.67	0.14	-9.44e-03	-0.10	1.13	0.46
Affective_touch	0.49	-0.13	0.29	0.13	1.97	0.59
Blood_Sugar	0.29	0.21	0.18	-0.23	3.51	0.75
Burp	-0.03	0.73	-0.03	-1.21e-03	1.00	0.50
Cough	7.56e-03	0.69	7.00e-03	-0.01	1.00	0.52
Sneeze	0.07	0.54	7.09e-03	0.20	1.30	0.55
Wind	0.08	0.46	0.07	0.16	1.37	0.62
Hungry	-0.16	0.04	0.65	0.03	1.13	0.60
Thirsty	-4.84e-03	0.01	0.49	0.11	1.11	0.70
Breathing	0.14	0.05	0.44	0.03	1.24	0.70
Heart	0.11	-0.07	0.40	0.08	1.32	0.78
Bruise	0.31	0.25	0.34	-0.25	3.74	0.59
Pain	0.29	0.01	0.32	0.20	2.67	0.61
Temp	0.19	0.13	0.29	0.16	2.91	0.67
Muscles	0.20	0.13	0.26	0.16	3.24	0.69
Sex_arousal	0.18	0.02	0.21	0.15	2.81	0.82
Defecate	0.10	0.11	2.77e-03	0.66	1.10	0.46
Urinate	-0.11	0.13	0.24	0.52	1.65	0.54
Taste	0.23	0.06	0.17	0.30	2.60	0.69
Vomit	0.21	0.18	0.11	0.25	3.26	0.71

The 4 latent factors (oblimin rotation) accounted for 37.94% of the total variance of the original data (MR2 = 10.46%, MR3 = 10.17%, MR1 = 10.08%, MR4 = 7.23%).

Sample 1b

efa4 <- parameters::factor_analysis(data1b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR2	MR3	MR4	Complexity	Uniqueness
Urinate	0.69	-0.18	0.13	-0.04	1.22	0.54
Defecate	0.68	-0.13	0.04	0.06	1.10	0.55
Breathing	0.56	0.03	0.06	7.96e-03	1.03	0.63
Hungry	0.55	0.15	-7.15e-03	0.02	1.16	0.60
Thirsty	0.52	0.05	0.14	-0.01	1.16	0.62
Sex_arousal	0.50	0.14	-0.04	0.20	1.49	0.56
Temp	0.46	0.23	0.16	-0.01	1.77	0.52
Pain	0.46	0.43	-0.05	-0.03	2.03	0.49
Taste	0.40	0.21	-0.07	0.12	1.77	0.69
Heart	0.40	0.06	-0.04	0.15	1.36	0.76
Muscles	0.35	0.24	0.07	0.09	2.07	0.65
Vomit	0.32	0.06	0.23	0.11	2.20	0.68
Tickle	0.02	0.69	0.01	0.04	1.01	0.47
Itch	-0.02	0.67	0.18	0.03	1.16	0.41
Bruise	0.07	0.44	0.15	0.06	1.32	0.65
Blood_Sugar	-0.08	0.38	0.16	0.09	1.61	0.77
Affective_touch	0.38	0.38	-0.13	0.05	2.26	0.63
Sneeze	0.05	0.05	0.71	0.01	1.02	0.41
Cough	0.04	0.02	0.71	0.04	1.01	0.43
Wind	7.52e-03	-0.04	-5.18e-03	0.94	1.00	0.15
Burp	-0.01	0.27	0.17	0.47	1.91	0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Sample 2

efa4 <- parameters::factor_analysis(data2, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR4	MR2	Complexity	Uniqueness
Urinate	0.62	0.12	0.03	-0.16	1.21	0.55
Sex_arousal	0.58	0.02	0.03	-3.44e-04	1.01	0.64
Hungry	0.50	-0.07	0.02	0.25	1.52	0.65
Defecate	0.50	0.25	0.04	-0.12	1.62	0.56
Thirsty	0.47	-6.03e-03	0.07	0.14	1.22	0.70
Breathing	0.38	6.09e-03	-0.02	0.31	1.93	0.71
Temp	0.35	0.26	0.06	0.11	2.13	0.63
Vomit	0.33	0.22	0.06	0.05	1.84	0.72
Heart	0.32	-0.08	0.02	0.28	2.10	0.80
Taste	0.28	0.28	0.10	-0.06	2.32	0.73
Cough	-7.64e-03	0.69	0.04	0.09	1.04	0.46
Sneeze	0.15	0.58	0.02	-0.06	1.16	0.57
Burp	0.12	0.42	0.09	0.15	1.51	0.63
Wind	0.19	0.42	0.07	0.08	1.52	0.63
Tickle	0.02	-0.01	0.83	-0.06	1.01	0.35
Itch	-0.11	0.23	0.51	0.19	1.80	0.54
Affective_touch	0.24	-0.19	0.37	0.19	2.89	0.72
Bruise	-0.06	0.23	0.05	0.53	1.40	0.58
Blood_Sugar	-0.04	-7.08e-03	0.22	0.40	1.56	0.75
Pain	0.20	0.23	0.02	0.36	2.33	0.63
Muscles	0.19	0.21	-0.01	0.27	2.75	0.75

The 4 latent factors (oblimin rotation) accounted for 36.65% of the total variance of the original data (MR1 = 12.95%, MR3 = 10.42%, MR4 = 6.75%, MR2 = 6.52%).

Sample 3

efa4 <- parameters::factor_analysis(data3, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Cough	0.88	-0.08	-0.04	0.11	1.05	0.29
Sneeze	0.74	0.05	1.42e-03	-0.03	1.01	0.41
Vomit	0.55	0.12	0.04	-0.09	1.17	0.58
Burp	0.55	0.08	0.22	-0.13	1.47	0.48
Wind	0.53	0.10	0.21	-0.25	1.84	0.48
Temp	0.46	0.34	-0.06	0.09	1.97	0.52
Taste	0.29	0.18	0.10	0.12	2.37	0.74
Hungry	-8.85e-03	0.70	-0.01	-0.09	1.03	0.54
Urinate	0.08	0.66	-0.12	0.11	1.15	0.52
Thirsty	-0.08	0.64	-0.03	0.14	1.14	0.62
Defecate	0.10	0.63	-3.50e-05	-0.07	1.08	0.53
Breathing	-0.10	0.49	0.28	-0.07	1.76	0.66
Sex_arousal	0.17	0.44	0.17	-0.22	2.20	0.58
Pain	0.21	0.41	0.08	0.18	2.00	0.58
Heart	0.05	0.28	0.24	0.08	2.23	0.77
Muscles	0.24	0.27	0.20	0.19	3.66	0.61
Tickle	0.02	-0.05	0.84	-0.04	1.01	0.31
Itch	0.03	0.02	0.74	0.11	1.05	0.39
Affective_touch	0.08	0.23	0.36	0.12	2.08	0.67
Bruise	0.17	0.07	0.29	0.45	2.08	0.55
Blood_Sugar	4.11e-03	0.11	0.35	0.43	2.07	0.60

The 4 latent factors (oblimin rotation) accounted for 45.69% of the total variance of the original data (MR1 = 15.50%, MR3 = 15.46%, MR2 = 11.02%, MR4 = 3.71%).

sample 4

efa4 <- parameters::factor_analysis(data4, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR4	MR2	MR3	Complexity	Uniqueness
Burp	0.69	-0.09	-4.33e-03	0.04	1.04	0.57
Cough	0.66	0.03	0.04	-1.08e-03	1.01	0.52
Wind	0.64	1.45e-03	0.03	0.02	1.01	0.56
Sneeze	0.52	0.07	0.07	-3.65e-03	1.07	0.64
Vomit	0.50	0.05	0.02	0.10	1.10	0.66
Pain	0.40	0.28	0.04	0.12	2.00	0.56
Muscles	0.39	0.23	0.12	0.07	1.91	0.60
Sex_arousal	0.29	0.28	0.08	0.02	2.18	0.70
Bruise	0.24	-0.13	0.20	0.23	3.47	0.77
Hungry	-0.12	0.62	0.04	0.17	1.25	0.56
Urinate	0.20	0.60	-0.02	-0.09	1.27	0.55
Thirsty	-0.10	0.53	0.09	0.15	1.30	0.66
Defecate	0.26	0.48	-0.04	-0.01	1.58	0.61
Temp	0.27	0.40	0.07	0.04	1.84	0.60
Taste	0.20	0.36	0.15	0.04	1.99	0.65
Itch	5.98e-03	-0.07	0.82	0.01	1.01	0.35
Tickle	8.91e-04	0.06	0.73	-0.05	1.03	0.47
Affective_touch	0.04	0.24	0.45	0.04	1.58	0.63
Blood_Sugar	0.06	-0.13	0.33	0.19	2.05	0.81
Breathing	2.50e-03	0.09	-0.03	0.78	1.03	0.34
Heart	0.06	-0.08	0.03	0.70	1.05	0.50

The 4 latent factors (oblimin rotation) accounted for 41.38% of the total variance of the original data (MR1 = 14.47%, MR4 = 10.83%, MR2 = 8.95%, MR3 = 7.14%).

Sample 5

efa4 <- parameters::factor_analysis(data5, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR4	MR2	MR3	Complexity	Uniqueness
Sex_arousal	0.53	0.02	-0.05	-0.10	1.09	0.73
Wind	0.52	0.14	0.06	0.06	1.19	0.60
Heart	0.50	-0.08	0.04	0.08	1.12	0.75
Taste	0.44	-0.03	-0.02	0.17	1.29	0.76
Muscles	0.42	-0.03	0.10	0.15	1.38	0.75
Hungry	0.38	0.06	-0.17	-3.75e-03	1.42	0.84
Burp	0.38	0.06	0.14	-0.17	1.78	0.80
Affective_touch	0.36	-4.99e-03	-0.07	0.03	1.09	0.87
Defecate	0.35	0.23	-0.14	0.03	2.08	0.76
Thirsty	0.34	0.22	-0.03	-0.10	1.98	0.78
Blood_Sugar	0.27	-0.18	0.05	-0.01	1.85	0.93
Sneeze	0.27	0.17	0.26	-0.08	2.91	0.75
Urinate	-0.05	0.72	-0.11	0.06	1.07	0.51
Cough	0.04	0.63	0.19	0.01	1.19	0.49
Pain	0.12	0.57	0.10	-0.10	1.21	0.59
Temp	0.29	0.30	-0.06	0.09	2.27	0.72
Itch	-0.04	0.01	0.94	0.03	1.01	0.12
Bruise	0.20	0.13	0.43	-0.03	1.64	0.68
Tickle	0.27	-0.18	0.38	0.04	2.30	0.75
Breathing	3.71e-03	3.31e-03	0.02	0.96	1.00	0.07
Vomit	0.25	0.26	-4.07e-03	0.29	2.97	0.67

The 4 latent factors (oblimin rotation) accounted for 33.72% of the total variance of the original data (MR1 = 12.07%, MR4 = 8.86%, MR2 = 7.16%, MR3 = 5.62%).

Sample 6

efa4 <- parameters::factor_analysis(data6, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR4	MR2	MR3	Complexity	Uniqueness
Burp	0.79	0.03	-0.05	0.02	1.01	0.38
Cough	0.66	0.06	0.17	-0.28	1.53	0.32
Wind	0.62	0.10	0.03	0.03	1.06	0.49
Bruise	0.58	0.04	-0.13	0.16	1.25	0.66
Breathing	0.55	-0.05	0.09	0.22	1.38	0.59
Muscles	0.49	0.20	0.06	0.09	1.42	0.53
Sneeze	0.49	0.08	0.32	-0.26	2.37	0.37
Blood_Sugar	0.49	0.02	-0.15	0.18	1.49	0.76
Heart	0.43	-0.04	7.61e-03	0.24	1.59	0.75
Vomit	0.43	0.09	0.16	-0.15	1.64	0.63
Temp	0.34	0.20	0.27	-0.09	2.72	0.52
Taste	0.25	0.20	0.18	0.10	3.13	0.69
Tickle	-0.05	0.89	-0.06	0.02	1.02	0.32
Itch	0.02	0.74	0.08	-0.10	1.06	0.37
Affective_touch	0.14	0.48	0.01	0.14	1.36	0.62
Pain	0.24	0.43	0.07	0.12	1.81	0.54
Sex_arousal	0.24	0.29	0.14	0.14	2.93	0.62
Defecate	0.04	0.02	0.79	5.46e-03	1.01	0.32
Urinate	-0.06	0.05	0.75	0.06	1.03	0.44
Hungry	0.13	0.12	0.28	0.47	1.94	0.54
Thirsty	0.10	0.03	0.38	0.44	2.09	0.58

The 4 latent factors (oblimin rotation) accounted for 47.47% of the total variance of the original data (MR1 = 20.28%, MR4 = 12.37%, MR2 = 10.55%, MR3 = 4.27%).

Sample 7a

efa4 <- parameters::factor_analysis(data7a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Itch	0.63	-0.12	0.22	0.07	1.34	0.46
Pain	0.59	-0.01	0.03	0.09	1.05	0.59
Tickle	0.52	-0.12	0.14	0.14	1.41	0.61
Muscles	0.49	0.13	-0.03	0.07	1.19	0.67
Affective_touch	0.45	0.13	0.08	-2.36e-03	1.24	0.70
Breathing	0.41	0.26	-0.16	-0.03	2.03	0.76
Bruise	0.38	0.35	-0.04	-0.05	2.05	0.68
Temp	0.33	-0.04	0.17	0.20	2.22	0.71
Blood_Sugar	0.33	0.33	0.02	-0.13	2.30	0.75
Heart	0.30	0.30	-0.06	-0.22	2.88	0.83
Burp	-0.04	0.66	0.20	0.07	1.22	0.39
Wind	-0.01	0.61	0.08	0.15	1.17	0.50
Cough	0.03	0.11	0.76	-7.60e-04	1.05	0.31
Sneeze	0.07	5.91e-03	0.70	8.40e-05	1.02	0.47
Vomit	0.11	0.25	0.25	0.18	3.18	0.67
Urinate	3.91e-04	0.06	2.65e-04	0.57	1.02	0.65
Defecate	0.06	0.07	0.04	0.56	1.06	0.61
Hungry	0.30	5.52e-03	-0.23	0.36	2.69	0.72
Sex_arousal	0.13	0.23	1.59e-03	0.34	2.09	0.71
Thirsty	0.24	0.12	-0.10	0.27	2.69	0.79
Taste	0.17	0.12	0.19	0.21	3.52	0.75

The 4 latent factors (oblimin rotation) accounted for 36.52% of the total variance of the original data (MR1 = 13.12%, MR3 = 8.41%, MR2 = 7.80%, MR4 = 7.18%).

Sample 7b

efa4 <- parameters::factor_analysis(data7b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Sneeze	0.68	0.07	0.04	-0.16	1.14	0.53
Cough	0.67	0.06	0.04	-0.07	1.04	0.51
Burp	0.65	-0.07	-0.02	0.19	1.21	0.50
Wind	0.53	-0.01	2.01e-03	0.20	1.28	0.59
Vomit	0.48	0.05	0.05	-7.54e-03	1.05	0.72
Temp	0.25	0.23	0.15	0.09	2.91	0.67
Muscles	0.25	0.24	0.06	0.20	3.05	0.65
Sex_arousal	0.24	0.08	0.21	0.16	3.04	0.70
Tickle	-0.06	0.75	-0.02	-3.98e-03	1.01	0.50
Itch	0.05	0.72	0.02	-0.02	1.01	0.43
Affective_touch	0.07	0.42	0.01	0.14	1.31	0.70
Pain	0.13	0.34	0.09	0.21	2.15	0.63
Urinate	-0.04	5.49e-03	0.75	-6.11e-03	1.01	0.46
Defecate	0.06	-1.14e-03	0.69	-0.04	1.02	0.50
Thirsty	-5.07e-03	-0.04	0.35	0.24	1.83	0.79
Taste	0.16	0.13	0.26	0.17	3.07	0.69
Breathing	0.02	0.11	0.09	0.47	1.19	0.68
Heart	0.03	0.12	-0.01	0.45	1.16	0.74
Blood_Sugar	0.08	0.05	0.02	0.38	1.14	0.79
Bruise	0.23	0.10	8.46e-03	0.34	1.98	0.70
Hungry	-3.08e-03	0.02	0.28	0.32	1.97	0.75

The 4 latent factors (oblimin rotation) accounted for 36.96% of the total variance of the original data (MR1 = 12.41%, MR3 = 9.28%, MR2 = 8.10%, MR4 = 7.16%).

Sample 7c

efa4 <- parameters::factor_analysis(data7c, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR3	MR1	MR4	MR2	Complexity	Uniqueness
Cough	0.73	-0.04	0.04	0.04	1.02	0.44
Sneeze	0.73	-4.69e-03	-0.03	9.49e-03	1.00	0.48
Wind	0.56	0.10	0.09	-2.95e-03	1.12	0.58
Vomit	0.50	0.04	-0.01	-6.16e-03	1.02	0.74
Burp	0.47	0.13	0.12	0.08	1.34	0.60
Urinate	0.10	0.56	-0.05	4.96e-03	1.09	0.65
Defecate	0.24	0.55	-0.11	-0.03	1.47	0.60
Hungry	-0.16	0.54	0.14	0.03	1.31	0.69
Sex_arousal	0.09	0.45	0.08	0.04	1.15	0.71
Thirsty	-0.06	0.39	0.20	0.03	1.56	0.76
Temp	0.18	0.27	0.21	0.05	2.82	0.71
Taste	0.12	0.25	0.22	0.05	2.53	0.77
Bruise	0.11	-0.09	0.54	0.12	1.26	0.61
Muscles	0.13	0.06	0.46	-0.10	1.32	0.72
Breathing	0.04	0.22	0.42	-0.08	1.62	0.71
Blood_Sugar	0.02	-0.12	0.41	0.10	1.30	0.81
Pain	0.02	0.29	0.39	0.06	1.90	0.64
Heart	0.08	0.07	0.31	0.05	1.29	0.84
Tickle	-0.01	0.06	-0.06	0.85	1.02	0.30
Itch	0.05	-0.08	0.07	0.78	1.04	0.34
Affective_touch	-0.05	0.29	0.12	0.31	2.37	0.75

The 4 latent factors (oblimin rotation) accounted for 35.87% of the total variance of the original data (MR3 = 11.03%, MR1 = 9.41%, MR4 = 7.85%, MR2 = 7.57%).

Sample 8a

efa4 <- parameters::factor_analysis(data8a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Sneeze	0.69	-0.03	-0.01	0.03	1.01	0.53
Cough	0.68	-0.07	0.03	0.04	1.04	0.54
Burp	0.55	0.18	-2.31e-04	-0.09	1.26	0.61
Wind	0.48	0.08	0.07	-0.03	1.09	0.70
Vomit	0.42	0.09	-0.04	0.13	1.31	0.74
Taste	0.20	0.19	0.18	0.04	3.08	0.79
Bruise	-0.03	0.65	-9.02e-03	0.04	1.01	0.58
Affective_touch	0.07	0.55	-0.03	-0.02	1.05	0.67
Blood_Sugar	0.15	0.45	-0.09	0.09	1.40	0.70
Itch	0.07	0.41	0.14	0.04	1.30	0.71
Muscles	0.02	0.39	0.15	0.02	1.29	0.77
Thirsty	5.49e-03	-0.05	0.67	-0.02	1.01	0.58
Hungry	0.06	-0.10	0.49	0.17	1.35	0.67
Pain	-0.04	0.31	0.41	0.02	1.89	0.67
Urinate	0.13	0.12	0.39	-2.22e-03	1.43	0.73
Temp	0.03	0.16	0.36	0.09	1.58	0.75
Defecate	0.24	0.09	0.34	-0.04	2.00	0.72
Heart	-0.02	0.04	-0.04	0.65	1.01	0.58
Breathing	0.08	-6.25e-03	0.12	0.49	1.17	0.66
Sex_arousal	0.06	0.19	0.17	0.22	3.04	0.78

The 4 latent factors (oblimin rotation) accounted for 32.61% of the total variance of the original data (MR1 = 10.53%, MR3 = 9.06%, MR2 = 8.24%, MR4 = 4.78%).

Sample 8b

efa4 <- parameters::factor_analysis(data8b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR4	MR2	MR3	Complexity	Uniqueness
Burp	0.65	0.07	6.06e-03	-0.10	1.07	0.56
Cough	0.63	-0.06	0.03	0.11	1.09	0.56
Sneeze	0.57	0.02	-8.72e-03	0.13	1.10	0.61
Wind	0.54	0.02	-0.03	1.35e-04	1.01	0.71
Defecate	0.34	0.11	0.28	-0.13	2.48	0.70
Vomit	0.27	0.18	0.06	0.13	2.39	0.77
Bruise	-0.06	0.64	-0.04	0.01	1.02	0.63
Affective_touch	0.08	0.55	-0.07	-0.08	1.11	0.67
Blood_Sugar	0.09	0.50	-0.03	0.02	1.07	0.71
Muscles	-2.84e-04	0.41	0.12	0.05	1.19	0.79
Itch	0.15	0.40	0.02	0.13	1.52	0.69
Pain	0.05	0.29	0.26	0.02	2.07	0.78
Taste	0.09	0.24	0.21	0.11	2.76	0.79
Thirsty	0.03	-0.10	0.71	-0.04	1.05	0.51
Hungry	-4.08e-03	0.04	0.51	0.17	1.23	0.65
Temp	-0.09	0.22	0.39	0.10	1.87	0.77
Urinate	0.19	0.18	0.33	-0.11	2.49	0.73
Breathing	0.06	-0.03	-0.01	0.70	1.02	0.50
Heart	-0.01	0.07	0.04	0.50	1.06	0.72
Sex_arousal	0.06	0.17	0.17	0.30	2.35	0.76

The 4 latent factors (oblimin rotation) accounted for 32.02% of the total variance of the original data (MR1 = 9.96%, MR4 = 9.20%, MR2 = 7.36%, MR3 = 5.49%).

Sample 9

efa4 <- parameters::factor_analysis(data9, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Burp	0.78	-0.01	-9.54e-03	-0.06	1.01	0.45
Wind	0.71	-0.03	-7.51e-03	0.03	1.01	0.51
Cough	0.69	0.02	0.01	0.04	1.01	0.47
Sneeze	0.63	-1.87e-03	0.02	0.11	1.07	0.51
Vomit	0.42	0.06	0.06	0.19	1.49	0.63
Temp	0.30	0.24	0.03	0.20	2.75	0.60
Sex_arousal	0.27	0.21	0.03	0.18	2.76	0.68
Taste	0.24	0.19	0.12	0.13	3.07	0.72
Breathing	0.04	0.62	-0.03	0.02	1.02	0.58
Hungry	-0.09	0.59	-0.01	0.15	1.17	0.60
Thirsty	-0.06	0.51	-0.01	0.23	1.42	0.60
Heart	0.08	0.47	0.06	-0.08	1.16	0.74
Pain	0.18	0.41	0.09	0.06	1.52	0.62
Muscles	0.34	0.37	0.06	7.80e-04	2.05	0.57
Blood_Sugar	0.08	0.31	0.28	-0.16	2.68	0.76
Bruise	0.24	0.29	0.28	-0.19	3.70	0.66
Tickle	-0.04	-0.04	0.79	0.05	1.02	0.41
Itch	0.01	-0.01	0.77	-0.02	1.00	0.41
Affective_touch	0.05	0.19	0.34	0.11	1.85	0.71
Urinate	0.02	0.09	0.03	0.69	1.04	0.41
Defecate	0.11	0.02	0.05	0.67	1.06	0.43

The 4 latent factors (oblimin rotation) accounted for 42.44% of the total variance of the original data (MR1 = 14.92%, MR3 = 11.60%, MR2 = 8.42%, MR4 = 7.50%).

Sample 10a

efa4 <- parameters::factor_analysis(data10a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR4	MR2	Complexity	Uniqueness
Hungry	0.62	-0.04	0.15	0.04	1.13	0.51
Sex_arousal	0.58	0.06	0.01	0.03	1.03	0.60
Pain	0.55	0.04	0.09	-8.09e-03	1.06	0.61
Thirsty	0.55	0.09	0.07	0.03	1.10	0.58
Temp	0.36	0.11	0.19	0.09	1.85	0.65
Affective_touch	0.28	-0.05	0.13	0.26	2.51	0.75
Defecate	0.06	0.67	-7.30e-03	0.02	1.02	0.49
Sneeze	-0.18	0.53	0.29	0.06	1.85	0.59
Vomit	0.13	0.53	-0.05	0.09	1.20	0.61
Urinate	0.40	0.47	-0.12	-0.02	2.09	0.52
Burp	2.70e-03	0.46	0.17	0.08	1.34	0.64
Breathing	0.02	0.02	0.71	0.01	1.00	0.45
Heart	0.17	7.44e-03	0.59	0.01	1.16	0.51
Muscles	0.15	0.07	0.47	0.07	1.31	0.60
Itch	-5.06e-03	-0.02	-0.02	0.78	1.00	0.43
Tickle	-0.01	0.02	-0.01	0.73	1.00	0.47
Bruise	8.63e-03	0.04	0.06	0.38	1.08	0.81
Wind	0.10	0.08	0.10	0.26	1.82	0.82

The 4 latent factors (oblimin rotation) accounted for 40.83% of the total variance of the original data (MR1 = 12.34%, MR3 = 10.09%, MR4 = 9.29%, MR2 = 9.12%).

Sample 10b

efa4 <- parameters::factor_analysis(data10b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR4	MR3	MR1	MR2	Complexity	Uniqueness
Breathing	0.55	-0.14	0.10	0.14	1.32	0.62
Muscles	0.55	0.16	0.02	-0.02	1.18	0.59
Heart	0.48	-0.01	0.04	0.08	1.07	0.71
Taste	0.42	0.20	-0.05	0.09	1.55	0.67
Temp	0.38	0.05	0.38	-0.06	2.09	0.56
Pain	0.36	0.07	0.32	-0.04	2.07	0.63
Cough	0.30	0.27	-0.07	0.16	2.62	0.69
Defecate	0.07	0.60	0.07	-0.06	1.07	0.58
Wind	-0.03	0.57	0.05	0.17	1.19	0.57
Burp	0.12	0.39	-0.08	0.27	2.13	0.64
Vomit	0.03	0.31	0.18	0.11	1.90	0.77
Sneeze	0.19	0.31	-0.06	0.23	2.64	0.70
Hungry	4.71e-03	0.04	0.63	0.07	1.03	0.55
Thirsty	0.21	-0.02	0.50	0.11	1.44	0.55
Urinate	-0.03	0.35	0.42	-0.03	1.96	0.62
Affective_touch	0.08	0.07	0.39	0.18	1.59	0.69
Sex_arousal	0.17	0.32	0.34	-0.10	2.64	0.61
Itch	0.15	0.09	-0.10	0.54	1.30	0.59
Bruise	-0.07	-0.03	0.23	0.53	1.39	0.67
Tickle	0.06	0.15	0.06	0.48	1.26	0.63
Blood_Sugar	0.05	-0.12	0.20	0.32	2.07	0.84

The 4 latent factors (oblimin rotation) accounted for 35.81% of the total variance of the original data (MR4 = 10.18%, MR3 = 9.19%, MR1 = 9.17%, MR2 = 7.27%).

Sample 10c

efa4 <- parameters::factor_analysis(data10c, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR2	MR4	MR3	Complexity	Uniqueness
Hungry	0.67	-0.11	-0.15	0.20	1.34	0.53
Thirsty	0.63	0.08	-0.04	-0.05	1.06	0.59
Sex_arousal	0.57	0.14	0.14	-0.07	1.28	0.53
Pain	0.54	0.07	0.14	0.07	1.20	0.57
Urinate	0.49	0.04	0.35	-0.08	1.89	0.52
Temp	0.35	0.04	0.35	0.15	2.37	0.56
Muscles	0.33	0.10	0.33	-0.09	2.35	0.66
Affective_touch	0.28	0.18	-0.02	0.13	2.15	0.81
Itch	-0.05	0.79	-0.01	-0.06	1.02	0.41
Bruise	0.07	0.55	-0.07	0.18	1.29	0.63
Taste	0.28	0.42	-0.14	0.09	2.10	0.67
Tickle	0.07	0.40	0.28	-0.06	1.90	0.66
Breathing	0.20	0.25	0.22	0.09	3.16	0.72
Wind	-0.05	-0.04	0.69	0.09	1.05	0.53
Sneeze	0.29	0.02	0.35	0.23	2.72	0.57
Defecate	0.25	0.20	0.28	0.04	2.85	0.68
Heart	0.04	0.12	0.22	0.04	1.66	0.90
Vomit	0.15	-0.01	0.17	0.04	2.11	0.92
Burp	0.03	-0.03	0.03	0.75	1.01	0.41
Blood_Sugar	-0.25	0.20	0.17	0.35	2.97	0.82
Cough	0.05	0.26	0.07	0.34	2.00	0.75

The 4 latent factors (oblimin rotation) accounted for 35.93% of the total variance of the original data (MR1 = 13.46%, MR2 = 8.82%, MR4 = 8.11%, MR3 = 5.54%).

Sample 11a

efa4 <- parameters::factor_analysis(data11a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR2	MR4	MR3	Complexity	Uniqueness
Cough	0.78	-0.02	-1.87e-03	0.15	1.07	0.34
Heart	0.57	-0.10	-0.04	-0.22	1.36	0.69
Breathing	0.56	0.14	0.07	-0.22	1.50	0.59
Sneeze	0.44	0.16	0.07	0.19	1.73	0.64
Wind	0.40	0.40	0.02	0.14	2.25	0.54
Temp	0.35	-0.08	0.32	-0.08	2.21	0.76
Vomit	0.32	0.11	0.17	-0.05	1.83	0.80
Defecate	8.51e-03	0.76	0.08	-0.09	1.05	0.38
Urinate	-0.08	0.67	0.19	-5.85e-03	1.19	0.45
Burp	0.35	0.45	-0.16	0.27	2.91	0.53
Sex_arousal	0.37	0.37	0.10	-0.10	2.28	0.60
Muscles	0.27	0.28	0.11	0.16	2.97	0.70
Thirsty	-0.13	0.16	0.74	0.08	1.19	0.36
Taste	0.14	0.12	0.49	-0.03	1.31	0.64
Blood_Sugar	0.30	-0.19	0.47	6.77e-03	2.09	0.68
Hungry	0.06	0.13	0.45	-0.07	1.25	0.72
Affective_touch	0.31	-0.23	0.36	0.08	2.82	0.74
Pain	0.11	0.29	0.33	0.12	2.51	0.65
Itch	-0.04	-0.04	-5.51e-03	0.72	1.01	0.49
Tickle	0.05	-0.03	0.12	0.52	1.13	0.69
Bruise	0.06	-0.17	0.18	0.41	1.79	0.78

The 4 latent factors (oblimin rotation) accounted for 39.22% of the total variance of the original data (MR1 = 12.83%, MR2 = 10.45%, MR4 = 9.59%, MR3 = 6.35%).

Sample 11b

efa4 <- parameters::factor_analysis(data11b, n = 4, rotation = "oblimin", sort = TRUE)

Warning in GPFoblq(A, Tmat = Tmat, normalize = normalize, eps = eps, maxit =
maxit, : convergence not obtained in GPFoblq. 1000 iterations used.

plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR2	MR1	MR4	MR3	Complexity	Uniqueness
Hungry	0.75	-0.09	0.03	0.06	1.05	0.46
Thirsty	0.74	0.07	-0.19	0.06	1.16	0.43
Breathing	0.71	0.08	-0.19	0.14	1.26	0.42
Urinate	0.67	0.01	0.16	-0.21	1.31	0.47
Pain	0.65	6.54e-03	0.35	0.10	1.61	0.28
Defecate	0.51	0.18	0.20	-0.28	2.22	0.47
Muscles	0.41	0.20	0.04	0.31	2.40	0.52
Temp	0.38	0.38	0.19	-0.02	2.45	0.43
Heart	0.37	0.12	0.19	0.06	1.79	0.70
Burp	0.08	0.83	-0.20	0.20	1.27	0.26
Sneeze	0.05	0.69	0.21	-0.32	1.64	0.26
Cough	0.02	0.66	0.09	-0.01	1.04	0.49
Wind	0.11	0.61	-0.02	0.06	1.10	0.55
Itch	-0.22	0.59	0.35	0.10	2.01	0.45
Vomit	0.12	0.51	0.22	-0.13	1.61	0.50
Tickle	0.10	0.26	0.50	0.05	1.62	0.50
Taste	0.14	0.19	0.47	-0.07	1.60	0.57
Affective_touch	0.10	0.12	0.41	0.35	2.29	0.60
Sex_arousal	0.22	0.17	0.41	0.14	2.24	0.57
Blood_Sugar	0.05	0.02	0.16	0.55	1.18	0.67
Bruise	0.27	0.17	0.11	0.50	1.92	0.48

The 4 latent factors (oblimin rotation) accounted for 51.95% of the total variance of the original data (MR2 = 19.08%, MR1 = 17.78%, MR4 = 9.34%, MR3 = 5.75%).

CFA

We discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant.

Ambiguous: - Temp - Vomit - Affective_touch - Sex_arousal - Taste

Lower-Level Factors

All data

m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"

m4 <- "
Sustenance =~ Hungry + Thirsty + Urinate + Defecate + Muscles + Pain
Skin =~ Itch + Bruise
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m5 <- "
Sustenance =~ Hungry + Thirsty + Pain
UrinateDefecate =~ Urinate + Defecate
Skin =~ Itch + Bruise + Muscles
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m6 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"

m7 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
CoughSneeze =~ Cough + Sneeze
WindBurp =~ Wind + Burp
"

is_converged <- function(m) {
  ev <- eigen(lavaan::lavTech(m, "cov.lv")[[1]], symmetric=TRUE, only.values=TRUE)$values
  if (any(ev < -1 * .Machine$double.eps^(3/4))) {
    "Not converged"
  } else {
    ""
  }
}

compare_mods <- function(data) {
  
  fit_model <- function(f, data) {
    if(!"Cough" %in% names(data)) {
      f <- str_remove(f, fixed("Cough +"))
    }
    lavaan::cfa(f, data = data)
  }
  
  
  cfa1 <- fit_model(m1, data = data)
  cfa4 <- fit_model(m4, data = data)
  cfa5 <- fit_model(m5, data = data)
  cfa6 <- fit_model(m6, data = data)
  cfa7 <- fit_model(m7, data = data)
  
  anova(cfa1, cfa4, cfa5, cfa6, cfa7) |> 
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(cfa1, cfa4, cfa5, cfa6, cfa7), is_converged))
}

display(compare_mods(data_all))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-171638	-171228	2195.20	56
cfa6	-167039	-166680	6805.73	62	< .001
cfa5	-162764	-162446	11091.21	67	< .001
cfa4	-157684	-157400	16178.96	71	< .001
cfa1	-144334	-144100	29540.73	77	< .001

m7_all <- lavaan::cfa(m7, data = data_all)
parameters::parameters(m7_all, standardize=TRUE) |>
  display()

# Loading

Link	Coefficient	SE	95% CI	z	p
HungryThirsty =~ Hungry	0.68	4.69e-03	(0.67, 0.69)	145.39	< .001
HungryThirsty =~ Thirsty	0.71	4.68e-03	(0.70, 0.72)	151.38	< .001
UrinateDefecate =~ Urinate	0.77	3.88e-03	(0.76, 0.78)	198.08	< .001
UrinateDefecate =~ Defecate	0.77	3.88e-03	(0.76, 0.78)	199.00	< .001
ItchBruise =~ Itch	0.50	5.65e-03	(0.49, 0.52)	89.31	< .001
ItchBruise =~ Bruise	0.67	5.89e-03	(0.66, 0.68)	113.21	< .001
MusclesPain =~ Muscles	0.71	4.10e-03	(0.70, 0.72)	173.39	< .001
MusclesPain =~ Pain	0.68	4.17e-03	(0.67, 0.69)	163.56	< .001
HeartBreathing =~ Breathing	0.79	5.12e-03	(0.78, 0.80)	153.73	< .001
HeartBreathing =~ Heart	0.59	5.06e-03	(0.58, 0.60)	116.02	< .001
CoughSneeze =~ Cough	0.80	3.47e-03	(0.80, 0.81)	230.97	< .001
CoughSneeze =~ Sneeze	0.75	3.62e-03	(0.74, 0.75)	206.40	< .001
WindBurp =~ Wind	0.77	3.55e-03	(0.76, 0.77)	216.49	< .001
WindBurp =~ Burp	0.82	3.42e-03	(0.81, 0.82)	239.00	< .001

# Correlation

Link	Coefficient	SE	95% CI	z	p
HungryThirsty ~~ UrinateDefecate	0.65	6.15e-03	(0.64, 0.66)	105.82	< .001
HungryThirsty ~~ ItchBruise	0.45	8.89e-03	(0.43, 0.47)	50.88	< .001
HungryThirsty ~~ MusclesPain	0.63	6.91e-03	(0.62, 0.65)	91.58	< .001
HungryThirsty ~~ HeartBreathing	0.64	6.90e-03	(0.63, 0.66)	93.27	< .001
HungryThirsty ~~ CoughSneeze	0.46	7.03e-03	(0.45, 0.47)	65.42	< .001
HungryThirsty ~~ WindBurp	0.43	7.02e-03	(0.42, 0.45)	61.80	< .001
UrinateDefecate ~~ ItchBruise	0.44	8.27e-03	(0.42, 0.46)	53.09	< .001
UrinateDefecate ~~ MusclesPain	0.65	6.14e-03	(0.64, 0.66)	105.65	< .001
UrinateDefecate ~~ HeartBreathing	0.53	6.71e-03	(0.52, 0.55)	79.61	< .001
UrinateDefecate ~~ CoughSneeze	0.56	5.95e-03	(0.55, 0.57)	94.13	< .001
UrinateDefecate ~~ WindBurp	0.51	6.10e-03	(0.50, 0.52)	83.40	< .001
ItchBruise ~~ MusclesPain	0.80	7.86e-03	(0.78, 0.81)	101.69	< .001
ItchBruise ~~ HeartBreathing	0.51	8.62e-03	(0.49, 0.53)	59.31	< .001
ItchBruise ~~ CoughSneeze	0.64	7.60e-03	(0.63, 0.65)	84.23	< .001
ItchBruise ~~ WindBurp	0.64	7.46e-03	(0.63, 0.66)	86.06	< .001
MusclesPain ~~ HeartBreathing	0.64	6.88e-03	(0.63, 0.65)	93.18	< .001
MusclesPain ~~ CoughSneeze	0.67	6.00e-03	(0.65, 0.68)	110.95	< .001
MusclesPain ~~ WindBurp	0.63	6.04e-03	(0.62, 0.65)	105.01	< .001
HeartBreathing ~~ CoughSneeze	0.52	6.73e-03	(0.50, 0.53)	76.91	< .001
HeartBreathing ~~ WindBurp	0.46	6.86e-03	(0.44, 0.47)	66.38	< .001
CoughSneeze ~~ WindBurp	0.72	4.80e-03	(0.71, 0.73)	150.93	< .001

Sample 1a

display(compare_mods(data1a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-1951	-1749	88.22	56
cfa6	-1925	-1748	126.29	62	< .001
cfa5	-1860	-1703	201.38	67	< .001	Not converged
cfa4	-1834	-1694	234.82	71	< .001
cfa1	-1661	-1546	420.24	77	< .001

Sample 1b

display(compare_mods(data1b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-1898	-1706	66.50	56
cfa6	-1823	-1654	154.19	62	< .001
cfa5	-1801	-1652	185.98	67	< .001
cfa4	-1755	-1621	240.22	71	< .001
cfa1	-1632	-1523	374.42	77	< .001

Sample 2

display(compare_mods(data2))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-4588	-4358	123.77	56
cfa6	-4565	-4363	158.82	62	< .001
cfa5	-4490	-4311	243.72	67	< .001
cfa4	-4426	-4266	315.85	71	< .001
cfa1	-4259	-4127	495.03	77	< .001

Sample 3

display(compare_mods(data3))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-3588	-3380	114.03	56
cfa6	-3508	-3325	206.45	62	< .001
cfa5	-3453	-3291	271.41	67	< .001
cfa4	-3419	-3275	313.05	71	< .001
cfa1	-3143	-3024	600.82	77	< .001

Sample 4

display(compare_mods(data4))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-3397	-3167	81.35	56
cfa6	-3321	-3120	168.76	62	< .001
cfa5	-3142	-2964	357.91	67	< .001
cfa4	-3069	-2910	438.71	71	< .001
cfa1	-2627	-2496	892.62	77	< .001

Sample 5

display(compare_mods(data5))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-96	49	101.12	56
cfa6	-101	26	107.73	62	0.359
cfa5	-102	10	116.58	67	0.115	Not converged
cfa4	-104	-3	123.39	71	0.146	Not converged
cfa1	-84	-1	154.66	77	< .001	Not converged

Sample 6

display(compare_mods(data6))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-4359	-4142	152.88	56
cfa6	-4146	-3956	376.94	62	< .001	Not converged
cfa5	-4061	-3893	472.34	67	< .001	Not converged
cfa4	-3865	-3714	676.52	71	< .001	Not converged
cfa1	-3699	-3575	854.37	77	< .001	Not converged

Sample 7a

display(compare_mods(data7a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-3268	-3059	106.90	56
cfa6	-3133	-2950	254.03	62	< .001	Not converged
cfa5	-3115	-2953	282.30	67	< .001	Not converged
cfa4	-3053	-2909	351.61	71	< .001	Not converged
cfa1	-2842	-2723	574.94	77	< .001	Not converged

Sample 7b

display(compare_mods(data7b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-11892	-11617	188.21	56
cfa6	-11703	-11462	388.89	62	< .001	Not converged
cfa5	-11583	-11370	518.70	67	< .001
cfa4	-11256	-11066	853.27	71	< .001	Not converged
cfa1	-10705	-10549	1416.33	77	< .001	Not converged

Sample 7c

display(compare_mods(data7c))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-1513	-1295	122.56	56
cfa6	-1429	-1237	219.24	62	< .001
cfa5	-1359	-1190	298.55	67	< .001
cfa4	-1296	-1144	370.09	71	< .001
cfa1	-1060	-935	617.75	77	< .001

Sample 8a

display(compare_mods(data8a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-5137	-4889	115.89	56
cfa6	-5103	-4885	161.81	62	< .001
cfa5	-5041	-4849	233.76	67	< .001
cfa4	-4930	-4758	352.62	71	< .001
cfa1	-4591	-4449	703.81	77	< .001

Sample 8b

display(compare_mods(data8b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-1738	-1531	83.61	56
cfa6	-1711	-1530	122.14	62	< .001
cfa5	-1683	-1523	160.21	67	< .001
cfa4	-1630	-1487	221.58	71	< .001
cfa1	-1486	-1368	377.91	77	< .001

Sample 9

display(compare_mods(data9))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-137262	-136870	2003.31	56
cfa6	-133150	-132807	6126.47	62	< .001
cfa5	-129698	-129394	9589.14	67	< .001
cfa4	-125505	-125233	13790.36	71	< .001
cfa1	-114370	-114146	24936.87	77	< .001

Sample 10a

display(compare_mods(data10a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-1824	-1632	101.11	45
cfa6	-1826	-1654	109.52	50	0.135
cfa5	-1822	-1671	123.48	55	0.016
cfa4	-1752	-1619	200.78	59	< .001
cfa1	-1665	-1556	300.09	65	< .001

Sample 10b

display(compare_mods(data10b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-3965	-3734	185.78	56
cfa6	-3959	-3755	204.27	62	0.005
cfa5	-3938	-3759	234.70	67	< .001	Not converged
cfa4	-3875	-3714	306.23	71	< .001
cfa1	-3710	-3578	482.94	77	< .001

Sample 10c

display(compare_mods(data10c))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-432	-303	67.23	56
cfa6	-438	-324	73.36	62	0.409
cfa5	-443	-343	78.32	67	0.421	Not converged
cfa4	-448	-358	81.66	71	0.502	Not converged
cfa1	-434	-360	107.51	77	< .001	Not converged

Sample 11a

display(compare_mods(data11a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-219	-88	62.63	56
cfa6	-208	-93	86.39	62	< .001
cfa5	-179	-78	124.90	67	< .001	Not converged
cfa4	-171	-80	141.41	71	0.002
cfa1	-99	-24	224.81	77	< .001

Sample 11b

display(compare_mods(data11b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
cfa7	-612	-471	133.75	56
cfa6	-597	-474	159.97	62	< .001	Not converged
cfa5	-589	-480	178.10	67	0.003	Not converged
cfa4	-587	-489	188.60	71	0.033	Not converged
cfa1	-506	-425	281.32	77	< .001	Not converged

Conclusion: The 7th factor solution was preferred in most datasets, including by indices that penalize increased number of parameters (such as BIC).

Higher Order Factors

All data

m7h1 <- paste0(m7, "
              Interoception =~ HungryThirsty + UrinateDefecate + ItchBruise + MusclesPain + HeartBreathing + CoughSneeze + WindBurp
              ")

m7h2 <- paste0(m7, "
              Expulsion =~ CoughSneeze + WindBurp
              Interoception =~ HungryThirsty + UrinateDefecate + HeartBreathing
              ")

m7h3 <- paste0(m7, "
              Valenced =~ MusclesPain + HeartBreathing + ItchBruise
              Expulsion =~ CoughSneeze + WindBurp
              Homeostasis =~ HungryThirsty + UrinateDefecate
              ")

compare_mods_h <- function(data) {
  
  fit_model_h <- function(f, data) {
    if(!"Cough" %in% names(data)) {
      f <- str_remove(f, fixed("Cough +"))
    }
    lavaan::cfa(f, data = data)
  }
  
  m7 <- fit_model_h(m7, data = data)
  m7h2 <- fit_model_h(m7h2, data = data)
  m7h1 <- fit_model_h(m7h1, data = data)
  m7h3 <- fit_model_h(m7h3, data = data)
  
  anova(m7, m7h1, m7h2, m7h3) |>
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(m7, m7h1, m7h2, m7h3), is_converged))
}

display(compare_mods_h(data_all))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-171638	-171228	2195.20	56
m7h2	-170853	-170527	2999.77	66	< .001
m7h3	-169568	-169250	4286.71	67	< .001
m7h1	-167097	-166805	6763.53	70	< .001

No evidence for higher order factors.

Sample 1a

display(compare_mods_h(data1a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1951	-1749	88.22	56
m7h2	-1930	-1770	128.74	66	< .001
m7h3	-1904	-1747	157.51	67	< .001
m7h1	-1873	-1729	193.71	70	< .001

Sample 1b

display(compare_mods_h(data1b))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1898	-1706	66.50	56
m7h2	-1901	-1748	83.62	66	0.072
m7h3	-1855	-1706	131.76	67	< .001	Not converged
m7h1	-1841	-1704	151.59	70	< .001

Sample 2

display(compare_mods_h(data2))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-4588	-4358	123.77	56
m7h2	-4573	-4390	158.33	66	< .001
m7h3	-4563	-4385	170.60	67	< .001
m7h1	-4536	-4372	203.47	70	< .001

Sample 3

display(compare_mods_h(data3))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-3588	-3380	114.03	56
m7h2	-3585	-3420	136.81	66	0.012
m7h3	-3581	-3419	143.30	67	0.011
m7h1	-3497	-3348	233.30	70	< .001

Sample 4

display(compare_mods_h(data4))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-3397	-3167	81.35	56
m7h2	-3344	-3161	154.44	66	< .001
m7h3	-3335	-3157	165.03	67	0.001
m7h1	-3282	-3118	223.71	70	< .001

Sample 5

display(compare_mods_h(data5))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-96	49	101.12	56		Not converged
m7h2	-95	21	122.05	66	0.022	Not converged
m7h3	-95	18	124.10	67	0.152	Not converged
m7h1	-100	3	124.72	70	0.893	Not converged

Sample 6

display(compare_mods_h(data6))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-4359	-4142	152.88	56		Not converged
m7h2	-4329	-4156	202.72	66	< .001	Not converged
m7h3	-4329	-4161	204.07	67	0.246	Not converged
m7h1	-4311	-4157	227.94	70	< .001	Not converged

Sample 7a

display(compare_mods_h(data7a))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-3268	-3059	106.90	56		Not converged
m7h2	-3260	-3094	134.89	66	0.002	Not converged
m7h3	-3255	-3093	141.68	67	0.009	Not converged
m7h1	-3208	-3059	195.06	70	< .001	Not converged

Sample 7b

display(compare_mods_h(data7b))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.
Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-11892	-11617	188.21	56		Not converged
m7h2	-11876	-11658	223.53	66	< .001	Not converged
m7h3	-11872	-11659	229.97	67	0.011	Not converged
m7h1	-11770	-11574	337.40	70	< .001	Not converged

Sample 7c

display(compare_mods_h(data7c))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1513	-1295	122.56	56
m7h2	-1499	-1325	156.56	66	< .001
m7h3	-1483	-1313	175.42	67	< .001
m7h1	-1437	-1281	226.58	70	< .001

Sample 8a

display(compare_mods_h(data8a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-5137	-4889	115.89	56
m7h2	-5120	-4923	152.67	66	< .001
m7h3	-5103	-4911	171.44	67	< .001
m7h1	-5007	-4830	273.77	70	< .001

Sample 8b

display(compare_mods_h(data8b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1738	-1531	83.61	56
m7h2	-1729	-1565	112.49	66	0.001
m7h3	-1727	-1567	116.28	67	0.052
m7h1	-1710	-1562	139.61	70	< .001

Sample 9

display(compare_mods_h(data9))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-137262	-136870	2003.31	56
m7h2	-136603	-136292	2681.71	66	< .001
m7h3	-135349	-135045	3938.35	67	< .001
m7h1	-133209	-132929	6083.88	70	< .001

Sample 10a

display(compare_mods_h(data10a))

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1824	-1632	101.11	45
m7h2	-1798	-1648	146.92	55	< .001
m7h3	-1794	-1647	153.35	56	0.011	Not converged
m7h1	-1778	-1645	174.82	59	< .001	Not converged

Sample 10b

display(compare_mods_h(data10b))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-3965	-3734	185.78	56
m7h2	-3938	-3753	233.42	66	< .001
m7h3	-3912	-3732	260.95	67	< .001
m7h1	-3899	-3734	279.76	70	< .001

Sample 10c

display(compare_mods_h(data10c))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-432	-303	67.23	56		Not converged
m7h2	-444	-341	75.22	66	0.630	Not converged
m7h3	-442	-341	79.80	67	0.032	Not converged
m7h1	-444	-351	83.57	70	0.287	Not converged

Sample 11a

display(compare_mods_h(data11a))

Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated ov variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-219	-88	62.63	56
m7h2	-220	-115	82.41	66	0.031
m7h3	-213	-111	91.40	67	0.003	Not converged
m7h1	-211	-117	99.05	70	0.054

Sample 11b

display(compare_mods_h(data11b))

Warning: lavaan->lav_object_post_check():  
   covariance matrix of latent variables is not positive definite ; use 
   lavInspect(fit, "cov.lv") to investigate.

Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative
Warning: lavaan->lav_object_post_check():  
   some estimated lv variances are negative

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-612	-471	133.75	56		Not converged
m7h2	-598	-486	167.45	66	< .001
m7h3	-580	-471	187.06	67	< .001	Not converged
m7h1	-555	-454	218.25	70	< .001	Not converged

Model Performance

rbind(
  performance::performance(m7_all, metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "All"),
  performance::performance(update(m7_all, data = data1a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1a"),
  performance::performance(update(m7_all, data = data1b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1b"),
  performance::performance(update(m7_all, data = data2), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m7_all, data = data3), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m7_all, data = data4), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m7_all, data = data5), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m7_all, data = data6), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(update(m7_all, data = data7a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7a"),
  performance::performance(update(m7_all, data = data7b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7b"),
  performance::performance(update(m7_all, data = data7c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7c"),
  performance::performance(update(m7_all, data = data8a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8a"),
  performance::performance(update(m7_all, data = data8b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8b"),
  performance::performance(update(m7_all, data = data9), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 9"),
  performance::performance(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data10a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10a"),
  performance::performance(update(m7_all, data = data10b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10b"),
  performance::performance(update(m7_all, data = data10c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10c"),
  performance::performance(update(m7_all, data = data11a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 11a"),
  performance::performance(update(m7_all, data = data11b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 11b")
) |>
  display()

Chi2	RMSEA	CFI	SRMR	Sample
2195.20	0.03	0.98	0.02	All
88.22	0.04	0.98	0.03	Sample 1a
66.50	0.02	0.99	0.03	Sample 1b
123.77	0.04	0.97	0.03	Sample 2
114.03	0.04	0.98	0.03	Sample 3
81.35	0.02	0.99	0.02	Sample 4
101.12	0.08	0.86	0.07	Sample 5
152.88	0.05	0.97	0.03	Sample 6
106.90	0.04	0.97	0.04	Sample 7a
188.21	0.03	0.98	0.02	Sample 7b
122.56	0.04	0.97	0.03	Sample 7c
115.89	0.03	0.98	0.02	Sample 8a
83.61	0.03	0.98	0.03	Sample 8b
2003.31	0.04	0.98	0.02	Sample 9
82.59	0.05	0.97	0.03	Sample 10a
185.78	0.05	0.95	0.03	Sample 10b
67.23	0.04	0.96	0.06	Sample 10c
62.63	0.03	0.98	0.06	Sample 11a
133.75	0.10	0.90	0.08	Sample 11b

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

# Add empirical variables
add_empirical <- function(data, df) {
  x <- data |>
    mutate(
      HungryThirsty = (df$Hungry + df$Thirsty) / 2,
      UrinateDefecate = (df$Urinate + df$Defecate) / 2,
      MusclesPain = (df$Muscles + df$Pain) / 2,
      HeartBreathing = (df$Heart + df$Breathing) / 2,
      ItchBruise = (df$Itch + df$Bruise) / 2,
      WindBurp = (df$Wind + df$Burp) / 2)
  if("Cough" %in% names(df)) {
    mutate(x, CoughSneeze = (df$Cough + df$Sneeze) / 2)
  } else {
    mutate(x, CoughSneeze = df$Sneeze)
  }
}

# Deal with missing data
temp <- predict(m7_all, newdata=data7c)
scores7c <- data.frame(matrix(ncol = 7, nrow = nrow(ega_scores_7c))) 
scores7c[!is.na(ega_scores_7c$EGA1), ] <- temp
names(scores7c) <- names(as.data.frame(temp))

scores <- list(
  sample1a = cbind(ega_scores_1a, data_addprefix(as.data.frame(predict(m7_all, newdata=data1a)), "CFA_")) |>
    add_empirical(data1a) |> 
    mutate(Sample = "Sample 1a"),
  sample1b = cbind(ega_scores_1b, data_addprefix(as.data.frame(predict(m7_all, newdata=data1b)), "CFA_")) |>
    add_empirical(data1b) |> 
    mutate(Sample = "Sample 1b"),
  sample2 = cbind(ega_scores_2, data_addprefix(as.data.frame(predict(m7_all, newdata=data2)), "CFA_")) |>
    add_empirical(data2) |> 
    mutate(Sample = "Sample 2"),
  sample3 = cbind(ega_scores_3, data_addprefix(as.data.frame(predict(m7_all, newdata=data3)), "CFA_")) |>
    add_empirical(data3) |> 
    mutate(Sample = "Sample 3"),
  sample4 = cbind(ega_scores_4, data_addprefix(as.data.frame(predict(m7_all, newdata=data4)), "CFA_")) |>
    add_empirical(data4) |> 
    mutate(Sample = "Sample 4"),
  sample5 = cbind(ega_scores_5, data_addprefix(as.data.frame(predict(m7_all, newdata=data5)), "CFA_")) |>
    add_empirical(data5) |> 
    mutate(Sample = "Sample 5"),
  sample6 = cbind(ega_scores_6, data_addprefix(as.data.frame(predict(m7_all, newdata=data6)), "CFA_")) |>
    add_empirical(data6) |> 
    mutate(Sample = "Sample 6"),
  sample7a = cbind(ega_scores_7a, data_addprefix(as.data.frame(predict(m7_all, newdata=data7a)), "CFA_")) |>
    add_empirical(data7a) |> 
    mutate(Sample = "Sample 7a"),
  sample7b = cbind(ega_scores_7b, data_addprefix(as.data.frame(predict(m7_all, newdata=data7b)), "CFA_")) |>
    add_empirical(data7b) |> 
    mutate(Sample = "Sample 7b"),
  sample7c = cbind(ega_scores_7c, data_addprefix(scores7c, "CFA_")) |> 
    add_empirical(data7c) |> 
    mutate(Sample = "Sample 7c"),
  sample8a = cbind(ega_scores_8a, data_addprefix(as.data.frame(predict(m7_all, newdata=data8a)), "CFA_")) |>
    add_empirical(data8a) |> 
    mutate(Sample = "Sample 8a"),
  sample8b = cbind(ega_scores_8b, data_addprefix(as.data.frame(predict(m7_all, newdata=data8b)), "CFA_")) |>
    add_empirical(data8b) |> 
    mutate(Sample = "Sample 8b"),
  sample9 = cbind(ega_scores_9, data_addprefix(as.data.frame(predict(m7_all, newdata=data9)), "CFA_")) |>
    add_empirical(data9) |> 
    mutate(Sample = "Sample 9"),
  sample10a = cbind(ega_scores_10a, data_addprefix(as.data.frame(predict(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data_all), newdata=data10a)), "CFA_")) |>
    add_empirical(data10a) |> 
    mutate(Sample = "Sample 10a", CFA_CoughSneeze = NA),
  sample10b = cbind(ega_scores_10b, data_addprefix(as.data.frame(predict(m7_all, newdata=data10b)), "CFA_")) |>
    add_empirical(data10b) |> 
    mutate(Sample = "Sample 10b"),
  sample10c = cbind(ega_scores_10c, data_addprefix(as.data.frame(predict(m7_all, newdata=data10c)), "CFA_")) |>
    add_empirical(data10c) |> 
    mutate(Sample = "Sample 10c"),
  sample11a = cbind(ega_scores_11a, data_addprefix(as.data.frame(predict(m7_all, newdata=data11a)), "CFA_")) |>
    add_empirical(data11a) |> 
    mutate(Sample = "Sample 11a"),
  sample11b = cbind(ega_scores_11b, data_addprefix(as.data.frame(predict(m7_all, newdata=data11b)), "CFA_")) |>
    add_empirical(data11b) |> 
    mutate(Sample = "Sample 11b")
)

save(scores, file = "../data/scores.RData")